package hos.map.tool;

import java.util.List;

import hos.map.coordinate.LatLng;
import hos.map.distance.BearingDistanceCache;
import hos.map.distance.LengthUnit;
import hos.map.geometry.Circle;

/**
 * <p>Title: CalculatorUtils </p>
 * <p>Description:  </p>
 * <p>Company: www.mapuni.com </p>
 *
 * @author : 蔡俊峰
 * @version : 1.0
 * @date : 2023-12-10 22:31
 */
public class CalculatorUtils {
    /**
     * 目标点是否包含在多边形内部
     */
    public static boolean isPolygonContains(List<LatLng> polygon, LatLng point) {
        if (polygon == null || polygon.size() < 3) {
            return false;
        }
        int size = polygon.size();
        double ax = 0;
        double ay = 0;
        double bx = polygon.get(size - 1).getLongitude() - point.getLongitude();
        double by = polygon.get(size - 1).getLatitude() - point.getLatitude();
        int depth = 0;

        for (int i = 0; i < polygon.size(); i++) {
            ax = bx;
            ay = by;
            bx = polygon.get(i).getLongitude() - point.getLongitude();
            by = polygon.get(i).getLatitude() - point.getLatitude();
            if (ay < 0 && by < 0) continue; // both "up" or both "down"
            if (ay > 0 && by > 0) continue; // both "up" or both "down"
            if (ax < 0 && bx < 0) continue; // both points on left
            double lx = ax - ay * (bx - ax) / (by - ay);
            if (lx == 0) return true; // point on edge
            if (lx > 0) depth++;
        }
        return (depth & 1) == 1;
    }

    /**
     * 返回以北极为顶点的三角形的带符号面积。根据Todhunter的《球面三角学》第71页第103节第2点，
     * 从“给定两条边和夹角的球面三角形的面积”推导出的公式。
     * The arguments named 'tan' are tan((pi/2 - x)/2).
     */
    private static double polarTriangleArea(double tan1, double lng1, double tan2, double lng2) {
        double deltaLng = lng1 - lng2;
        double t = tan1 * tan2;
        return 2 * Math.atan2(t * Math.sin(deltaLng), 1 + t * Math.cos(deltaLng));
    }

    /**
     * 计算集合的面积
     *
     * @param path 集合的路径
     * @return 面积
     */
    public static double area(List<LatLng> path) {
        if (path == null || path.size() < 3) {
            return 0;
        }
        int size = path.size();
        double total = 0;
        double pi = MathUtil.getPI();
        LatLng prev = path.get(size - 1);
        double prevTanLat = Math.tan((pi / 2 - prev.getLongitudeInRad()) / 2);
        double prevLng = prev.getLatitudeInRad();
        // For each edge, accumulate the signed area of the triangle formed by the North Pole
        // and that edge ('polar triangle').
        for (LatLng latLng : path) {
            double tanLat = Math.tan((pi / 2 - latLng.getLongitudeInRad()) / 2);
            double lng = latLng.getLatitudeInRad();
            total += polarTriangleArea(tanLat, lng, prevTanLat, prevLng);
            prevTanLat = tanLat;
            prevLng = lng;
        }
        return total * (MathUtil.getEarthRadius() * MathUtil.getEarthRadius());
    }

    /**
     * 根据 Haversine 计算两点之间的距离
     *
     * @param p1 第一个点
     * @param p2 第二个点位
     * @return 距离
     */
    public static double distanceHaversine(LatLng p1, LatLng p2) {
        if (p1 == null || p2 == null) {
            return 0;
        }
        final double sinDLat = Math.sin((p2.getLatitudeInRad() - p1.getLatitudeInRad()) / 2);
        final double sinDLng = Math.sin((p2.getLongitudeInRad() - p1.getLongitudeInRad()) / 2);
        // Sides
        final double a = sinDLat * sinDLat +
                sinDLng * sinDLng *
                        Math.cos(p1.getLatitudeInRad()) *
                        Math.cos(p2.getLatitudeInRad());
        final double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
        return MathUtil.getEarthRadius() * c;
    }

    /**
     * 根据 Haversine 计算点位的偏移
     *
     * @param from            点位
     * @param distanceInMeter 偏移的距离
     * @param bearing         偏移的角度
     * @return 偏移后的位置
     */
    public static LatLng offsetHaversine(LatLng from, double distanceInMeter, double bearing) {
        if (from == null || distanceInMeter == 0) {
            return null;
        }
        if (bearing < -180 || bearing > 180) {
            return null;
        }
        final double h = MathUtil.degToRadian(bearing);
        final double a = distanceInMeter / MathUtil.getEarthRadius();
        final double lat2 = Math.asin(Math.sin(from.getLatitudeInRad()) * Math.cos(a) +
                Math.cos(from.getLatitudeInRad()) * Math.sin(a) * Math.cos(h));
        final double lng2 = from.getLongitudeInRad() +
                Math.atan2(Math.sin(h) * Math.sin(a) * Math.cos(from.getLatitudeInRad()),
                        Math.cos(a) - Math.sin(from.getLatitudeInRad()) * Math.sin(lat2));
        return new LatLng(MathUtil.radianToDeg(lat2), MathUtil.radianToDeg(lng2), from.getType());
    }

    /**
     * 根据 Vincenty 计算两点之间的距离
     *
     * @param p1 第一个点
     * @param p2 第二个点位
     * @return 距离
     */
    public static BearingDistanceCache distanceVincenty(LatLng p1, LatLng p2) {
        if (p1 == null || p2 == null) {
            return null;
        }
        double lat1 = p1.getLatitude();
        double lat2 = p2.getLatitude();
        double lon1 = p1.getLongitude();
        double lon2 = p2.getLongitude();
        // Based on http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
        // using the "Inverse Formula" (section 4)

        int MAXITERS = 20;
        // Convert lat/long to radians
        lat1 *= Math.PI / 180.0;
        lat2 *= Math.PI / 180.0;
        lon1 *= Math.PI / 180.0;
        lon2 *= Math.PI / 180.0;

        double a = 6378137.0; // WGS84 major axis
        double b = 6356752.3142; // WGS84 semi-major axis
        double f = (a - b) / a;
        double aSqMinusBSqOverBSq = (a * a - b * b) / (b * b);

        double L = lon2 - lon1;
        double A = 0.0;
        double U1 = Math.atan((1.0 - f) * Math.tan(lat1));
        double U2 = Math.atan((1.0 - f) * Math.tan(lat2));

        double cosU1 = Math.cos(U1);
        double cosU2 = Math.cos(U2);
        double sinU1 = Math.sin(U1);
        double sinU2 = Math.sin(U2);
        double cosU1cosU2 = cosU1 * cosU2;
        double sinU1sinU2 = sinU1 * sinU2;

        double sigma = 0.0;
        double deltaSigma = 0.0;
        double cosSqAlpha = 0.0;
        double cos2SM = 0.0;
        double cosSigma = 0.0;
        double sinSigma = 0.0;
        double cosLambda = 0.0;
        double sinLambda = 0.0;

        double lambda = L; // initial guess
        for (int iter = 0; iter < MAXITERS; iter++) {
            double lambdaOrig = lambda;
            cosLambda = Math.cos(lambda);
            sinLambda = Math.sin(lambda);
            double t1 = cosU2 * sinLambda;
            double t2 = cosU1 * sinU2 - sinU1 * cosU2 * cosLambda;
            double sinSqSigma = t1 * t1 + t2 * t2; // (14)
            sinSigma = Math.sqrt(sinSqSigma);
            cosSigma = sinU1sinU2 + cosU1cosU2 * cosLambda; // (15)
            sigma = Math.atan2(sinSigma, cosSigma); // (16)
            double sinAlpha = (sinSigma == 0) ? 0.0 :
                    cosU1cosU2 * sinLambda / sinSigma; // (17)
            cosSqAlpha = 1.0 - sinAlpha * sinAlpha;
            cos2SM = (cosSqAlpha == 0) ? 0.0 :
                    cosSigma - 2.0 * sinU1sinU2 / cosSqAlpha; // (18)

            double uSquared = cosSqAlpha * aSqMinusBSqOverBSq; // defn
            A = 1 + (uSquared / 16384.0) * // (3)
                    (4096.0 + uSquared *
                            (-768 + uSquared * (320.0 - 175.0 * uSquared)));
            double B = (uSquared / 1024.0) * // (4)
                    (256.0 + uSquared *
                            (-128.0 + uSquared * (74.0 - 47.0 * uSquared)));
            double C = (f / 16.0) *
                    cosSqAlpha *
                    (4.0 + f * (4.0 - 3.0 * cosSqAlpha)); // (10)
            double cos2SMSq = cos2SM * cos2SM;
            deltaSigma = B * sinSigma * // (6)
                    (cos2SM + (B / 4.0) *
                            (cosSigma * (-1.0 + 2.0 * cos2SMSq) -
                                    (B / 6.0) * cos2SM *
                                            (-3.0 + 4.0 * sinSigma * sinSigma) *
                                            (-3.0 + 4.0 * cos2SMSq)));

            lambda = L +
                    (1.0 - C) * f * sinAlpha *
                            (sigma + C * sinSigma *
                                    (cos2SM + C * cosSigma *
                                            (-1.0 + 2.0 * cos2SM * cos2SM))); // (11)

            double delta = (lambda - lambdaOrig) / lambda;
            if (Math.abs(delta) < 1.0e-12) {
                break;
            }
        }
        float distance = (float) (b * A * (sigma - deltaSigma));
        BearingDistanceCache bearingDistanceCache = new BearingDistanceCache(distance);
        float initialBearing = (float) Math.atan2(cosU2 * sinLambda,
                cosU1 * sinU2 - sinU1 * cosU2 * cosLambda);
        initialBearing *= 180.0 / Math.PI;
        bearingDistanceCache.setInitialBearing(initialBearing);
        float finalBearing = (float) Math.atan2(cosU1 * sinLambda,
                -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda);
        finalBearing *= 180.0 / Math.PI;
        bearingDistanceCache.setFinalBearing(finalBearing);
        bearingDistanceCache.setLat1(lat1);
        bearingDistanceCache.setLat2(lat2);
        bearingDistanceCache.setLon1(lon1);
        bearingDistanceCache.setLon2(lon2);
        return bearingDistanceCache;
    }

    /**
     * 根据 Vincenty 计算点位的偏移
     *
     * @param from            点位
     * @param distanceInMeter 偏移的距离
     * @param bearing         偏移的角度
     * @return 偏移后的位置
     */
    public static LatLng offsetVincenty(LatLng from, double distanceInMeter, double bearing) {
        if (from == null || distanceInMeter == 0) {
            return null;
        }
        if (bearing < -180 || bearing > 180) {
            return null;
        }
        /// Equator radius in meter (WGS84 ellipsoid)
        double equatorRadius = 6378137.0;
        /// Polar radius in meter (WGS84 ellipsoid)
        double polarRadius = 6356752.314245;
        /// WGS84
        double flattening = 1 / 298.257223563;

        final double latitude = from.getLatitudeInRad();
        final double longitude = from.getLongitudeInRad();

        final double alpha1 = MathUtil.degToRadian(bearing);
        final double sinAlpha1 = Math.sin(alpha1);
        final double cosAlpha1 = Math.cos(alpha1);

        final double tanU1 = (1 - flattening) * Math.tan(latitude);
        final double cosU1 = 1 / Math.sqrt((1 + tanU1 * tanU1));
        final double sinU1 = tanU1 * cosU1;

        final double sigma1 = Math.atan2(tanU1, cosAlpha1);
        final double sinAlpha = cosU1 * sinAlpha1;
        final double cosSqAlpha = 1 - sinAlpha * sinAlpha;
        final double dfUSq = cosSqAlpha *
                (equatorRadius * equatorRadius - polarRadius * polarRadius) /
                (polarRadius * polarRadius);
        final double a = 1 +
                dfUSq / 16384 * (4096 + dfUSq * (-768 + dfUSq * (320 - 175 * dfUSq)));
        final double b = dfUSq / 1024 * (256 + dfUSq * (-128 + dfUSq * (74 - 47 * dfUSq)));

        double sigma = distanceInMeter / (polarRadius * a);
        double sigmaP = 2 * MathUtil.getPI();

        double sinSigma = 0.0;
        double cosSigma = 0.0;
        double cos2SigmaM = 0.0;
        double deltaSigma;
        double maxIterations = 200;

        do {
            cos2SigmaM = Math.cos(2 * sigma1 + sigma);
            sinSigma = Math.sin(sigma);
            cosSigma = Math.cos(sigma);
            deltaSigma = b *
                    sinSigma *
                    (cos2SigmaM +
                            b /
                                    4 *
                                    (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
                                            b /
                                                    6 *
                                                    cos2SigmaM *
                                                    (-3 + 4 * sinSigma * sinSigma) *
                                                    (-3 + 4 * cos2SigmaM * cos2SigmaM)));
            sigmaP = sigma;
            sigma = distanceInMeter / (polarRadius * a) + deltaSigma;
        } while (Math.abs(sigma - sigmaP) > 1e-12 && --maxIterations > 0);

        if (maxIterations == 0) {
            return null;
        }

        final double tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1;
        final double lat2 = Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,
                (1 - flattening) * Math.sqrt(sinAlpha * sinAlpha + tmp * tmp));

        final double lambda = Math.atan2(
                sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1);
        final double c =
                flattening / 16 * cosSqAlpha * (4 + flattening * (4 - 3 * cosSqAlpha));
        final double l = lambda -
                (1 - c) *
                        flattening *
                        sinAlpha *
                        (sigma +
                                c *
                                        sinSigma *
                                        (cos2SigmaM +
                                                c * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));

        double lon2 = longitude + l;
        // print("LA ${radianToDeg(lat2)}, LO ${radianToDeg(lon2)}");

        if (lon2 > MathUtil.getPI()) {
            lon2 = lon2 - 2 * MathUtil.getPI();
        }
        if (lon2 < -1 * MathUtil.getPI()) {
            lon2 = lon2 + 2 * MathUtil.getPI();
        }

        return new LatLng(MathUtil.radianToDeg(lat2), MathUtil.radianToDeg(lon2), from.getType());
    }


    /**
     * 计算给定路径在地球上的长度，以米为单位。
     */
    public static double distanceVincenty(List<LatLng> path) {
        if (path == null || path.size() < 2) {
            return 0;
        }
        double length = 0;
        int size = path.size();
        for (int i = 0; i < size - 1; i++) {
            LatLng current = path.get(i);
            LatLng next = path.get(i + 1);
            BearingDistanceCache bearingDistanceCache = distanceVincenty(current, next);
            if (bearingDistanceCache != null) {
                length += bearingDistanceCache.getDistance();
            }
        }
        return length;
    }

    /**
     * 计算给定路径在地球上的长度，
     */
    public static double distanceVincenty(LengthUnit unit, List<LatLng> path) {
        if (path == null || path.size() < 2) {
            return 0;
        }
        final double dist = distanceVincenty(path);
        // If the distance is NaN or infinite, return 0.0
        if (Double.isNaN(dist) || Double.isInfinite(dist)) {
            return 0;
        }
        return LengthUnit.Meter().to(unit, dist);
    }

    /**
     * 根据 Vincenty 计算两点之间的距离
     *
     * @param p1 第一个点
     * @param p2 第二个点位
     * @return 距离
     */
    public static BearingDistanceCache distanceVincenty(LengthUnit unit, LatLng p1, LatLng p2) {
        if (p1 == null || p2 == null) {
            return null;
        }
        BearingDistanceCache bearingDistanceCache = distanceVincenty(p1, p2);
        if (bearingDistanceCache == null) {
            return null;
        }
        final double dist = bearingDistanceCache.getDistance();
        // If the distance is NaN or infinite, return 0.0

        if (Double.isNaN(dist) || Double.isInfinite(dist)) {
            return null;
        }
        bearingDistanceCache.setDistance(LengthUnit.Meter().to(unit, dist));
        return bearingDistanceCache;
    }

    /**
     * 计算2个点位之间的角度
     *
     * @param p1 第一个点位
     * @param p1 第二个点位
     */
    public static double bearing(LatLng p1, LatLng p2) {
        if (p1 == null || p2 == null) {
            return 0;
        }
        final double diffLongitude = p2.getLongitudeInRad() - p1.getLongitudeInRad();

        final double y = Math.sin(diffLongitude);
        final double x = Math.cos(p1.getLatitudeInRad()) * Math.tan(p2.getLatitudeInRad()) -
                Math.sin(p1.getLatitudeInRad()) * Math.cos(diffLongitude);

        return MathUtil.radianToDeg(Math.atan2(y, x));
    }

    /**
     * 判断点位是否在圆内
     *
     * @param point 目标点位
     */
    public static boolean isPointInsideVincenty(final LatLng center, int radius, final LatLng point) {
        if (center == null || point == null) {
            return false;
        }
        BearingDistanceCache bearingDistanceCache = distanceVincenty(center, point);
        if (bearingDistanceCache == null) {
            return false;
        }
        final double dist = bearingDistanceCache.getDistance();
        return dist <= radius;
    }

    /**
     * 判断点位是否在圆内
     *
     * @param point 目标点位
     */
    public static boolean isPointInsideVincenty(final Circle center, final LatLng point) {
        if (center == null || point == null) {
            return false;
        }
        BearingDistanceCache bearingDistanceCache = distanceVincenty(center.getCenter(), point);
        if (bearingDistanceCache == null) {
            return false;
        }
        final double dist = bearingDistanceCache.getDistance();
        return dist <= center.getRadius();
    }

    /**
     * 计算给定路径在地球上的长度，以米为单位。
     */
    public static double distanceHaversine(List<LatLng> path) {
        if (path == null || path.size() < 2) {
            return 0;
        }
        double length = 0;
        int size = path.size();
        for (int i = 0; i < size - 1; i++) {
            LatLng current = path.get(i);
            LatLng next = path.get(i + 1);
            length += distanceHaversine(current, next);
        }
        return length;
    }

    /**
     * 计算给定路径在地球上的长度
     */
    public static double distanceHaversine(LengthUnit unit, List<LatLng> path) {
        if (path == null || path.size() < 2) {
            return 0;
        }
        final double dist = distanceHaversine(path);
        // If the distance is NaN or infinite, return 0.0
        if (Double.isNaN(dist) || Double.isInfinite(dist)) {
            return 0;
        }
        return LengthUnit.Meter().to(unit, dist);
    }

    /**
     * 计算两个点之间的距离
     */
    public static double distanceHaversine(LengthUnit unit, LatLng p1, LatLng p2) {
        if (p1 == null || p2 == null) {
            return 0;
        }
        final double dist = distanceHaversine(p1, p2);
        // If the distance is NaN or infinite, return 0.0

        if (Double.isNaN(dist) || Double.isInfinite(dist)) {
            return 0;
        }
        return LengthUnit.Meter().to(unit, dist);
    }

    /**
     * 判断点位是否在圆内
     *
     * @param point 目标点位
     */
    public static boolean isPointInsideHaversine(final LatLng center, int radius, final LatLng point) {
        if (center == null || point == null) {
            return false;
        }
        final double dist = distanceHaversine(center, point);
        return dist <= radius;
    }

    /**
     * 计算球面上点p和线段开始到结束之间的距离。
     */
    public static BearingDistanceCache distanceVincentyToLine(final LatLng p, final LatLng start, final LatLng end) {
        if (p == null || start == null || end == null) {
            return null;
        }
        if (start == end) {
            return distanceVincenty(end, p);
        }

        final double s0lat = p.getLatitudeInRad();
        final double s0lng = p.getLongitudeInRad();
        final double s1lat = start.getLatitudeInRad();
        final double s1lng = start.getLongitudeInRad();
        final double s2lat = end.getLatitudeInRad();
        final double s2lng = end.getLongitudeInRad();

        final double s2s1lat = s2lat - s1lat;
        final double s2s1lng = s2lng - s1lng;
        final double u = ((s0lat - s1lat) * s2s1lat + (s0lng - s1lng) * s2s1lng) / (s2s1lat * s2s1lat + s2s1lng * s2s1lng);
        if (u <= 0) {
            return distanceVincenty(p, start);
        }
        if (u >= 1) {
            return distanceVincenty(p, end);
        }
        final LatLng su = new LatLng(start.getLatitude() + u * (end.getLatitude() - start.getLatitude()),
                start.getLongitude() + u * (end.getLongitude() - start.getLongitude()), p.getType());
        return distanceVincenty(p, su);
    }

    /**
     * 计算球面上点p和线段开始到结束之间的距离。
     */
    public static double distanceHaversineToLine(final LatLng p, final LatLng start, final LatLng end) {
        if (p == null || start == null || end == null) {
            return 0;
        }
        if (start == end) {
            return distanceHaversine(end, p);
        }

        final double s0lat = p.getLatitudeInRad();
        final double s0lng = p.getLongitudeInRad();
        final double s1lat = start.getLatitudeInRad();
        final double s1lng = start.getLongitudeInRad();
        final double s2lat = end.getLatitudeInRad();
        final double s2lng = end.getLongitudeInRad();

        final double s2s1lat = s2lat - s1lat;
        final double s2s1lng = s2lng - s1lng;
        final double u = ((s0lat - s1lat) * s2s1lat + (s0lng - s1lng) * s2s1lng) / (s2s1lat * s2s1lat + s2s1lng * s2s1lng);
        if (u <= 0) {
            return distanceHaversine(p, start);
        }
        if (u >= 1) {
            return distanceHaversine(p, end);
        }
        final LatLng su = new LatLng(start.getLatitude() + u * (end.getLatitude() - start.getLatitude()),
                start.getLongitude() + u * (end.getLongitude() - start.getLongitude()), p.getType());
        return distanceHaversine(p, su);
    }
}
